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We investigate existence and stability of viscoelastic shock profiles for a class of 
planar models including the incompressible shear case studied by Antman and Malek- 
Madani. We establish that the resulting equations fall into the class of symmetrizable 
P I hyperbolic-parabolic systems, hence spectral stability implies linearized and nonlin- 

ear stability with sharp rates of decay. The new contributions are treatment of the 
• compressible case, formulation of a rigorous nonlinear stability theory, including verifi- 

cation of stability of small-amplitude Lax shocks, and the systematic incorporation in 
^ our investigations of numerical Evans function computations determining stability of 

^ large-amplitude and or nonclassical type shock profiles. 

^ 1 Introduction 

in 

In this paper, generalizing work of Antman and Malek-Madani |AMj in the incompressible 
shear flow case, we carry out the numerical and analytical study of the existence and 
stability of planar viscoelastic traveling waves in a 3d solid, for a simple prototypical elastic 
energy density, both for the general compressible and the incompressible shear flow case. 
<^ We establish that the resulting equations fall into the class of symmetrizable hyperbolic- 

^ parabolic systems studied in |MaZ21 lMaZ31 lMaZ4t \RZ\ IZ4] . hence spectral stability implies 

^ linearized and nonlinear stability with sharp rates of decay. This important point was 

^ previously left undecided, due to a lack of the necessary abstract stability framework. 

5-H The new contributions beyond what was done in |AM) are: treatment of the com- 

pressible case, consideration of large-amplitude waves (somewhat artificial given our simple 
choice of energy density; however, the methods used clearly generalize to more physically 
correct models), formulation of a rigorous nonlinear stability theory including verification 
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of stability of small-amplitude Lax waves, and the systematic incorporation in our investi- 
gations of numerical Evans function computations determining stability of large-amplitude 
and or nonclassical type shock profiles. For related analysis in various different settings, see 
IBHRZl iHLZl IHLyZj ICHNZl jBHZ, .BLZj. 



In preparation for future generalizations, we also discuss the case of phase-transitional 
viscoelasticity. It would be interesting to carry out similar analysis for more general classes 
of elastic energy density, as well as for the phase-transitional case which involves, at the 
technical level higher order dispersive terms relating to surface energy, and at the physical 
level, presumably, interesting new behaviors. 

Acknowledgment. Thanks to Stuart Antman and Constantine Dafermos for several 
helpful conversations, and to Stuart Antman for making available the working notes [A] . 



2 The equations of viscoelasticity 

The equations of isothermal viscoelasticity are given through the following balance of linear 
momentum: 

(2.1) - Vx ■ (dW{VO + Z{Vi, Vit)) = 0. 



Here, ^ : Q,x M_|_ — > denotes the deformation of a reference configuration C M'^ which 
models a viscoelastic body with constant temperature and density. A typical point in VL is 
denoted by X, so that the deformation gradient is given as: 

with the key constraint of: 

det F > 0. 



In (2.1) the operator Vjf stands for the divergence of an appropriate field. We use the 
convention that the divergence of a matrix field is taken row-wise. In what follows, we 
shall also use the matrix norm \F\ = (tr(F^F))^/^, which is induced by the inner product: 
Fi : F2 = tr(Ff F2). 

The mapping DW : M^^'^ — )■ M^^^ is the Piola-Kirchhoff stress tensor which, in agree- 
ment with the second law of thermodynamics, is expressed as the derivative of an elastic 
energy density W : M^^^ — > M_|_. The viscous stress tensor is given by the mapping 
Z : M^^^ X M^^^ — > M^^^, depending on the deformation gradient F and the velocity 
gradient Q = Ft = V£,t = Vv, where v = S^f 



The first order version of the inviscid part of ( |2.1[ ): 
(2.2) itt - Vx ■ (dW{VO) = 



is: 



(2.3) iF,v)t + Y,dx,{G,{F,v))=0. 



i=l 



2 



Above, {F, v) : Q 



1)12 



represents conserved quantities, while Gi : 



)12 



d12 



given by: 



-Gi{F, v) = v^ei e v^ei v^Ci 



d 



OF, 



ki 



■W{F) 



i = 1..3 



k=l 



are the fluxes, and e,- denotes the i-th coordinate vector of ^ 



2.1 The elastic energy density W 

The principle of material frame invariance imposes the following condition on W , with 
respect to the group SO (3) of proper rotations in M^: 



(2.4) 



W{RF) = W{F) VF G 



!)3x3 



Vi2 G S0(3). 



Also, the material consistency requires that: 



(2.5) 



W{F) +00 as det F ^ 0. 



In what follows, we shall restrict our attention to the class of isotropic materials, for whom 
the energy W satisfies additionally: 



(2.6) 



W{FR) = W{F) VF G 



p3x3 



yR G S0(3). 



Recall |Ba] that hyperbolicity of rt2.2^ is equivalent to rank-one convexity of W . 



(2.7) 



A particular example of W satisfying (2.4) and (2.6) is: 

1 



W^{F) 



1, 



F^F -Id|' 



F'F|^-2|F|^ + 3) 



and by a direct calculation, we obtain: 

DWo{F) = F{F^F-ld). 



Note that Wq in (2.7) is not quasiconvex (or polyconvex) as it is not globally rank-one 



convex. This follows by checking the Legendre-Hadamard condition. Indeed, for any A G 
1^3x3 Qj^g ^g^. Q2^^Wo{F) = {FA^F + FF'^A + AF^F - A) : A. Taking ^ = Id and 
F G skew we obtain a^^Wo(F) = |F|2 - 3 which is negative for |F| < ^/2>. 

On the other hand we see that (9^^Wo(^) = 2\sym{AR^)\'^ for R G S0(3). If rank A = 1 
then rank(Ai?'^) = 1 so syia{AR^) / 0. Therefore d\^Wo{R) > c\A\^ for every R G SO{2>) 
and every rank-one matrix A, with a uniform c > 0. This implies that Wq is rank-one 
convex in a neighborhood of 50(3). 

Notice also that Wq has quadratic growth close to 50(3). Indeed, write F = R + E, 



where for F close to 50(3) we have: R 



50(3) 



F and |F| = dist(F, 50(3)). Since F 



is orthogonal to the tangent space to 50(3) at R, we see that R E must be symmetric. 
Therefore: Wo{F) = WR'^E + E'^R + F'^Fp = l\2R^E + E'^ E\^ = |Fp + 0{\E\^). 



For other examples of W satisfying (2.4) and (2.6), see (3.7) and Appendix A.l 
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2.2 The viscous stress tensor Z 

The viscous stress tensor Z : M^^^ x M'^^^ — > M^^^ should be compatible with the following 
principles of continuum mechanics: balance of angular momentum, frame invariance, and 
the Claussius-Duhem inequality. That is, for every F,Q £ M'^^^ with det F 7^ 0, we require 
that: 



(i) skew ^Z{F, Q)) = 0, i.e. Z = FS with S symmetric. 

(ii) Z{RF, RtF + RQ) = RZ{F, Q) for every path of rotations : M+ — > S0{2,), 
i.e. in view of (i): S{RF, RKF + RQ) = S{R, Q) Vi? G SO{3) G skew. 



(iii) Z{F,Q) : Q > 0, i.e. in view of (i): S : sym(F^Q) > 0. 
Examples of Z satisfying the above are: 

Zi{F,Q) = 2FsYm{F^Q), 

Z2{F, Q) = 2(detF)sym(QF-i)F-i'^. 

We note that in the case of Z2, the related Cauchy stress tensor T2 = 2(deti^)^^^2-^"^ = 
2sym(QF"i) is the La grangian version of the stress tensor 2symVt' written in Eulerian co- 
ordinates. For incompressible fluids 2div(symVt;) = Af , giving the usual parabolic viscous 
regularization of the fluid dynamics evolutionary system. For more general viscous stress 



tensors, see Appendix A. 2 



2.3 An extension: the surface energy 



A phenomenological modification that is sometimes used is to replace (2.1) with 
(2.10) itt - Vx ■ [dW{VO + 2(Ve, V6) - f (V^O) = 0, 

where the surface energy £ is given by: 



Ji,fc:1...3 



for some convex density ^' : M^xSxs — y compatible with frame indifference (and 
isotropy). A typical example is ^o(G') = ||Gp, so that: 

(2.11) Soi^^O = Vx • V^e = AxF 

which is an extension of the Id case of [Sl]. 

Writing the variation of the energy J ^'(V^^) in the direction of a test function cj) £ 
C^{n,R^) we obtain: 



(2.12) / /^^(V^O : V^c/, = / {Vx-£(V^O)- 

Jn Jn 
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which justifies the last divergence term in (2.10). 

The addition of surface energy is motivated by the van der Waals/Cahn-HilUard ap- 
proach to the stationary equihbrium theory p H IZ81 ICGSl ISZ] . This would be an interesting 
direction for further investigation. 



2.4 Entropy 

After integrating (|2.1l) against on and then integrating by parts, we obtain: 



where we used that: 



^J{DW(yo + z{vc,vCt))n = o on on, 



a natural assumption following from either (Dirichlet) clamped boundary conditions (,^qq = 
const or else free (Neumann) conditions [DWCV^) + Z(V^, V^t)) = corresponding to 
the absence of stress on the boundary. 
Consequently: 

dt J 05*161' + Dw(yo^ = - j 2(yt V6) : V6 < 

by the Clausius-Duhem inequality, and we see that the integral / of the quantity: 



(2.13) 



7]iF,v) = -\v\'^ + WiF) 



along F = and v = ^t, is nonincreasing in time. In case of (2.10), using (2.12) with 
= we obtain that J fj is nonincreasing, for: 



Further, notice that t] : M}"^ 
Indeed, the scalar fields: 



defined in (2.13) is an entropy [D] associated to (2.3). 



d 



1 3 



dFik 



W{F) 



k=l 



qi{F,v) = -'. 

define the respective entropy fluxes, in the sense that: 

Vqi{F,v) = V7]{F,v)DGi{F,v). 



i = 1..3 
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3 The planar case 



We now restrict our attention to the interesting subclass of planar solutions, which are 
solutions in full 3d space that depend only on a single coordinate direction. Namely, we 
assume that the deformation ^ has the form: 

^{X) = X + Uiz), X = {x,y,z), U = {u,v,w) 

which yields the following structure of the deformation gradient: 







' 1 





Uz 




' 1 





ai 


(3.1) 


F = 





1 


Vz 







1 


0-2 












l+Wz _ 










as 



We shall denote V = {a,b) = (ai, 02, as, 6i, 62, ^s), where ai = n^,a2 = Vz,as = 1 + ^2 and 
61 = tif, 62 = Vt,b3 = Wt, with the constraint: 

(3.2) as > 0, 

corresponding to det F > in the region of physical feasibility of 

Writing W{a) = W{ 



1 ai 
1 aa 
as 



, we see that for all F as in (|3.1|) there holds: 
Vx • {DW{F)) = {DaW{a))z. 



That is, the planar equations inherit a vector-valued variational structure echoing the 
matrix-valued variational structure, and thus (|2.2|) has the following form: 



Vt + GiV)z = 
(3.3) G{V) = i-b, -DaW{a)f, DG{V) = 





-M 



-Ids 




M = DlW{a). 



It follows that strict hyperbolicity of (3.3) is equivalent to strict convexity of W with M 
having 3 distinct (positive) eigenvalues. Also, r]{V) = ^\b\'^ + W{a) is then a convex entropy: 



(3.4) 



V7]{V)DGiV) = Vq{V), q{V) = -b ■ DaWia). 



3.1 Energy density Wq and viscosities Zi 

By a straightforward calculation, we have: 



(3.5) 



Woia) 



\{\a\'-ir + l{al + al) 
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and: 

div(i:'Wo(F)) = (^Uzz + {UzlUzl"^ + 2UzWz)z, Vzz + {VzlUzl"^ + 2VzWz)z, 

(3.6) 2wzz + (1^7.1^ + Wz\Uz\^ + 2wl)zY 

= ((|a|V)., {\a\\2)z, ((|a|2-l)a3).)^ 
More generally, one may consider densities of the form: 

(3.7) WiF) = ^\F^F - Hp + C2(|F|2 - 3)^ + C3(|detF| - l)^. 

The C2 term contributes to DW{F) as: 4c2(|-Fp — 3)F in the planar case; this is Ac2{2wz + 
\Uz\'^)F, with divergence: 

4c2((u^|t/^|2 + 2uzWz)z, ivz\Uz\'^ + 2vzWz)z, {\Uz\'^ + 2wz + Wz\Uz\'^ + 2u;2)/ 



4c2(((|a|2 - l)ai)„ ((|a|2 - 1)02)., ((|a|2 - 1)03). 



T 



The C3 term contributes to DW{F) the term: 2(detF — l)cofF, for F with detF > 0. In 
the planar case, divergence of this term reads 2c3{0,0,Wzz)'^ = 2c3(0,0, (as)^)"^. 

Combining, we obtain the general form: 
(3.8) 

div(i:>W^(F)) = (^((^i|a|^ + ti2)ai)z, ((w|ap + fJ'2)a2)z, ((/^i|ap + /^2)a3 + (^3 - 1)03)2) , 
with m = 1 + 4c2, fi2 = — 4c2, fJ-3 = 2c3, corresponding to elastic potential 

W{a) = + ^Ai2|ap + ^/U3(a3 - 1)^ + C, 

where C is a constant. The above potential is strictly convex at the identity (a = (0,0, 1)) 
whenever /ii + ^2 > and 3/ii + ^2 + /£3> 0. It is a simple case of the general form 
W{a) = C7(|ap,a3) described in Appendix A.l Restricted to the incompressible planar case 
of Section 3.2.2 (3.8) recovers the class of equations studied in [AM| . 

Regarding the viscous tensors, we obtain: 

div(Zi(F,F)) = (uzzt + {Uz{2Wzt + {\Uz\^)t))z, Vzzt + {Vz{2Wzt + (1^7^^)*))^, 

(3.9) 2Wzzt + {\Uz?)zt + {Wz{2Wzt + {\Uz\^)t))z)' 

. T 

bi^z + 2aia ■ hz] , ( ^2,2 + 2a2a -hz] , ( 2030 • 6^ ) ^ 



, T 



(3.10) 



div(^2(F,i^)) = ( (^^) , (^^) , 

bi,z \ fb2,z\ 2(^) 
J \ / z^ \ / z 
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where a ■ bz = aibi z + 0262,2 + 0363,2. 



Hence, (2.1) with W and Z as in (2.7) and (2.9), has the hyperbolic-parabohc form: 
(3.11) Vt + G{V), = {B{V)Vz)z 



with Giy) as in (3.3) and: 

M = dIWq = diag(^|a|^ |a|^ |ap - 1^ + 2a (g) a 



in view of (3.5|). Further: 
(3.12) B : 





Bo^i 



So,i = diag(l,l,0) +2a(^a or ^0,2 = — diag( 1, 1, 2 



in case of Zi and Z2, respectively. Both tensors i?o,i are symmetric and positive definite on 
the entire physical region 03 > 0. 

3.2 The full and the restricted systems in hyperbolic— parabolic form 
3.2.1 Compressible viscoelasticity 



For the viscous stress tensor Zi, system (|3.11|) reads: 
(3.13) 



ai,t - bi,z 
0-2,1 - 62,2 
as,* - 63,2 



0, 
0, 
0, 



62, t 

63, * 



(|a|V) 
(|opa2) 



{bi^z + 2aia • bz 
(62,2 + 2a2a • bz 



{{\a\ - 1)03)2 = (2030 • 6^). 



For the viscous tensor Zo we have: 



(3.14) 



ai,t - 61,2 = 0, 
a2,t - 62,2 = 0, 
03,* - 63,2 = 0, 



62, * 

63, t 



{\a\'^ai)z = 
(|opa2)2 = 
((lap -1)03) 



61,2 

03 / z 

62,2 
03 



P3,2 

03 



3.2.2 The 2D incompressible shear case 



reduces to the following one (naturally, we now denote a = (ai, 02) and |ap = + alJT 



For an incompressible medium and a shear deformation where w = 0, the system (3.14) 

itu 

0, 



(3.15) 



ai,t - 61,2 
a2,t - 62,2 



0, 



bi,t - ((|ap + 1)01)2 = 61,22, 
b2,t - ((kP + 1)02)2 = 62,22, 



with an associated pressure of p = |ap whose gradient (0, 0, (|ap)2)"^ cancels the term 
— (|ap)2 in the 63 equation of (3.14). Note that the viscous stress tensor in this case reduces 
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to the Laplacian. Equations (3.15) are a special case of the equations studied in |AM] : they 
may be also recognized as the model for an elastic string. 
For the choice we obtain: 



(3.16) 



0.21 - b2,z 



0, 
0, 



6i,t- ((|a|2 + l)ai), 
b2,t-{(.\a\^ + I)a2)z 



{hi^z + 2aia • 
(62,2 + 2a2a • bz 



The incompressible model may be viewed as the formal limit as /i — t- +00 of a system with 
potential VFo(a) + /i3(a3 — 1)^, penalizing variations in density detF = 03. Operationally, 
this amounts to fixing 03 = 1 in a given (compressible) elastic potential and dropping the 
equation for 03, to obtain a reduced shear potential W{ai, 02) = Ty(ai, 02, 1) and equations 
whose first-order part have the same variational structure ( |3.3| ) as the full 3d system. 

3.2.3 The 2D compressible case 



Another reduced version of (3.14), restricted to the v — w plane is obtained by setting -u = 0. 



write |ap = 0^ + a|): 



This is an equally simple system as (3.15), but with essentially different structure (we now 



(3.17) 



a2,t - b2,z = 0, 
03,* - ^3,2 = 0, 



52,t 



(l«l'«2) 



bj^ 
as 



63,t- (dap -1)03). 



b3,z 



while for Zi, writing a-hz = 02^2,2 + 03^3,21 we have: 
(3.18) 

3.2.4 The ID cases 



a2,t - b2,z = 0, 
a3,t - bs^z = 0, 



b2,t - (|apa2)2 = (^2,2 + 2a2a • bz)^ , 
b3,t - ((|a|^ - 1)03)2 = (2a3a • bz)z ■ 



Taking v = w = 0, (3.14) further reduces to a model of the transverse unidirectional 



perturbations in a beam or string: 
(3.19) ai,t - bi,z = 0, 



bi,t - (a? + ai)z = bi^zz- 



Setting u = V = 0, ( |3.14[ ) yields the ID compressible model for longitudinal perturba- 

^3,2' 



tions in a viscoelastic rod: 
(3.20) a3,t - 63,2 = 



b3,t - (03 - as). 



«3 
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3.2.5 Extension: surface energy and higher-order dispersion 

We mention briefly the effects of modifying by the addition of surface energy term. In the 



planar, incompressible shear case, (|2.10|) with (|2.11|) becomes: 
(3.21) 

See [Slj for a corresponding treatment of the one-dimensional case. 



ai,t 

«2,t 



0, 
0, 



(ai + {of + al)ai) 
(a2 + {of + al)a2) 



bi,zz 

b2,zz 



Oil.zzz ) 
a2,zzz- 



3.3 Hyperbolic characteristics 
3.3.1 Compressible case 



Consid er th e inviscid version of (3.11): Vt + DG(y)Vz = 0. Using the block structure of 
DG in (3.3) we obtain that its eigenvalues are {±^mj}j^-^ with corresponding eigenvectors 
({rj, =F\/"^'"i)}j=i5 where nij (and rj) are the eigenvalues (and corresponding eigenvectors) 
of the symmetric matrix M. Also, since rrij are independent of 6, the linear degeneracy or 
genuine nonlinearity of the ziz^nij characteristic fields of DG is equivalent to the same 
properties of the mj characteristic fields of M. 

Using now the following formula, valid for 3x3 matrices: det{A + B) = detA + (cofj4) : 
B + (coiB) : A + detB, we obtain: 



(3.22) 



mi 

1712 



4|a| 



(2|a|2-l)2 + 8(a2 + a2) 
I (4|a|2 - 1 + ^(2|a|2-l)2 + 8(a? + ai)) 



Note that at a = (0, 0, 1) we have mi = m2 = 1 and ma = 2; hence DG is (nonstrictly) 
hyperbolic at Vb = (0, 0, 0, 6i, 62, 63). Further calculations show that, whenever defined: 

(i) ri = (02, — ai, 0)"^ and the eigenvalues it-^mi correspond to two linearly degenerate 
fields of DG. 



(iij r2 



-2aia3, — 2a2a3, 3|a| 



(iii) r3 = (-2aia3, - 20203, 3|a| 



2ai 
24 



1712) 



mzY' and in the vicinity of Vq the eigenvalues 



ib-y/m3 correspond to two genuinely nonlinear fields of DG. 
3.3.2 The 2D incompressible shear case 

bt-{DaW{a)), = b,,. 



The system (3.15) can be written as: 

at-bz = 0, 
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Its flux matrix depends on a - 
(3.23) «G,_, = 



(ai, 02) and has the form: 

, Mi_2 = DlWoia) 



|a|^ + l)Id2 + 2a<E)a, 



where Wo(a) = ||a|*^ + ^|a|2. We see that strict hyperbohcity, convexity of Wq and existence 
of strictly convex entropy are equivalently satisfied here. 

Calculating as before, DG1-2 has two genuinely nonlinear characteristic fields, with 
eigenvalues zby^l + 3\a'^\ and corresponding eigenvectors 

(ibai, ±02, ai\/l + 3|a|2, 02 ^/l + 3|ap)'^ 



in fast modes, and two linearly degenerate fields with eigenvalues zt-y/l + \a'^\ and eigenvec- 
tors (±02, ^ai, a2^/l + |a|2, — aiy^l + |a|2)"^ in slow modes. The linear degeneracy reflects 
the rotational degeneracy of the underlying system [FJ . 

3.3.3 The 2D compressible case 



The flux matrix in (3.17) depends on a = (02, as) and has the form: 



(3.24) 



DG2-3 





-M2-3 



-Ida 




M2-3 = diagi |ap, |ap - 1 ) + 2a 



and we flnd that DG2^3 tia-s two couples of eigenvalues {ziz^mj}j=2.3 with correspond- 



4|a|' 



a|2 - 1)2 -I- , 



ing eigenvectors (r^, =p ^nijrj), where: m2 = 2 

I (^4|a|2 - 1 + y^(2|a|2 - 1)2 + 8a^) , while = (-2a2a3, 3|a|2-2a|-mj)'^ (or ra = (1,0)"^ 

when a2 = 0). We see that the in the vicinity of (0, 1,62)^3) the matrix DG2-3 is strictly 
hyperbolic, the two eigenfields corresponding to zb-y/mjj are genuinely nonlinear. 

3.3.4 The ID incompressible case 



For system (3.19) the characteristic speeds are ±y'T+~3af, while for the system (3.20) they 
are ±-^303 — 1. Hence the second model is strictly hyperbolic for |a3| > l/\/3 and elliptic 
otherwise; this can be recognized as agreeing with certain phase-transitional viscoelasticity 
models, except that the region 03 < (where det -F < 0) is unphysical. 

4 Nonlinear stability framework 

We now briefly recall the general stability theory of \LA\ iRj IRZj , which reduces the question 



of nonlinear stability in (3.11 ) to veriflcation of an Evans function condition. Namely, given 



two endstates V- and Vj^ belonging to the regions of strict hyperbohcity of DG, we make 
a smooth change of coordinates V 1— >• SiV) with S : — > given by: S{V) = Dr](V) = 
DaW{a) © b. The system (3.11) is equivalent to: 



A''iS)St + A{S)S, = iBiS)S, 
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where with a shght abuse of notation we use S = S oV : [0, cxd) x 



V-^. Above: 



A = DG{V)A° 



-Ida 
-MQ 



B = B{V)A'^ = B{V), i° 



Q 
Ida 



where Q = Q(y) is defined as follows. In some open neighborhoods of and V+ (where 
M is positive definite) we set Q = M~^, in which case: 



M-i 
Ids 



In the region where M is negative definite, we set Q = Ida. In between the two above 
mentioned regions, Q is a smooth, symmetric and positive definite interpolation of the two 
matrix fields M and Ida • This construction allows us to treat also the case of profiles passing 
through elliptic regions, but with hyperbolic endstates (in a similar spirit as for the van der 
Waals gas dynamics examples mentioned in |MaZ41 IZ4j ) . 

We first check the validity of the structural conditions (Al)-(A3) of jiZ4] : 

(Al) AiV-) and ^(V+) are symmetric matrices. Aq is symmetric and positive definite 
(on the whole M^). Also, the 3x3 principal minor of A, corresponding to the purely 



hyperbolic part of the system (3.11), equals identically O3 hence it is always symmetric, as 
required. 

(A2) At the endstates V± there holds: no eigenvector of DG belongs to the kernel of B. 
In the region of strict hyperbolicity of DG this condition is equivalent to: no eigenvector of 
M is in the kernel of -Bo,i> readily satisfied. 

' Oa 



(A3) B has the required block structure B 



as in (3.12 



O3 

Oa Bo 

tion of the minor corresponding to the parabolic part of the system 
is uniformly elliptic in any region in V of the form: < 03 



The symmetriza- 



( |3.1ip : sym Bo^i = Bo,i 
< C in case of -Bq 2 > and 



> c(l + af 



in case of Bq^i (where c, C > are some uniform constants). 



We hence find that shock profiles of each of the planar systems considered in this paper 
satisfy conditions (Al)~(A3) of |Z4| defining the class of symmetrizable hyperbolic-parabolic 
systems and profiles to which the theory of nonlinear stability of viscous shock profiles 
developed in |MaZ2[ IMaZ3l IMaZ4l [H IH |RZ] apphes, provided: 

(i) the endstates V± lie in the region of strict hyperbolicity of DG, 

(ii) The profile {t^(-)} lies in some region where the chosen i?o,j is uniformly elliptic. 



We now validate the additional technical conditions (II0)-(II3) of [Z4j . Note that the 
remaining conditions (H4)-(H5) are needed only for the multi-dimensional systems, as they 
automatically hold for systems in 1 space dimension. 

(HO) G,B,Se C^. 

(HI) the shock speed s under consideration is non-zero (note that is the only eigenvalue 
of the 3x3 principal minor of DG, which indeed is Oa). As remarked in section [5| s 7^ 
for any profile with endstates belonging to the strict hyperbolicity region of DG. 
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(H2) s is distinct from the eigenvalues of DG{V± 



(H3) local to V{-), the set of traveling wave solutions to (3.11) connecting {V-,V+) 
(with thus determined speed s), forms a smooth finite-dimensional submanifold {t^^(-)} of 
Ci(M,M*^), parametrized hy 6 e B{0,r) C M^ and V° = V. 

4.1 The Evans condition 



Linearizing the hyperbolic-parabolic system (3.11 ) about its viscous shock solution of (3.11 ): 



V{z,t) = V{z-st), lim V{z) = V±, 

which satisfies: —sV+{G{V))z = {B{V)Vz)z, and further changing to co-moving coordinates 
z = z — st, we obtain the equivalent evolution equations: 

(4.1) Vt = CV := {BV,), - {GV),. 

Here Q and B are the following matrix fields depending on z: 

G{z) = DG{V{z)) - sId - DB{V{z)fV,{z), B{z) = B{V{z)), 

and converging asymptotically to values ^(±c«) = DG{V±) — sId and i3(±oo) = B{y±). 



Towards investigating stability of (4.1), one seeks eigenvalues A E C of £, that is solutions 



to the system CV = W written in its first-order form: 
(4.2) Z'{z,\) = A{z,\)Z{z,\). 

The augmented "phase variable" Z consists oi V = (a, h) and the derivative b' of its 
parabolic-like component. 

As shown in [G^ IZHl lMaZ3l IMaZ4[ IZl] . under conditions (HO), (HI), (H2) it is possible 



to define an analytic Evans function -D : {A G C;i?e A > 0} — > C associated with (4.2) 



and hence consequently associated with £ and with the original problem (3.11). We shall 



now briefly sketch this construction, for further details see e.g. |AGJt IGZl IZ4|, IHuZj . 



In the first step one observes that the complex matrix field A{z,\) G C"X" (4.2) is 



analytic in A and has an exponential decay to the respective .4-i-(A) as z — >• ±00 (uniformly 



in bounded A). The second step consists in proving that (4.2) on each of the half-lines 
(—00,0] and [0,00), is equivalent to: 

Z'{z) = A^{\)Z{z), z<0 and Z' {z) = A+{\) Z {z) , z > 0, 

under change of variables Z{z) = P-{z, \)Z{z) for z < 0, and Z{z) = P+{z, X)Z{z) for 
z >0. Existence of such (non-unique) analytic in A and invertible matrix fields P±{z, A) G 
£^NxN ^ decaying exponentially to Id as 2; — )• ±00, is achieved by a conjugation lemma |Z4| . 

Further, denote by {Z^{X)}i=i,,k the (analytic in A) basis of the stable space S of ^+(A), 
and likewise let {Z~ {X)}i=k+i..N be the basis of the unstable space U of ^_(A), where the 
consistency of the dimensions follows from assumptions (HI), (H2). Define: 

Z+iz,X) = P+{z,X)Z+iX), z>0 and Zr {z, X) = P-{z, X)Z- (X), z<0. 
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Clearly, given any Zq G spanjZj (zq, A)}j=i..fc, zq > 0, there exists a solution to (4.2) on 
[20,00) decaying exponentially to as z — )• 00, and with initial data Z{zo) = Zq. It has the 
property that Z {z, A ) G span{Z^ (z, X)}i=i..k for all z > zq- A similar assertion of backward 



resolvability of (4.2) is true for Zq S span{Z~ {zq, X)}i=k+i..N, -^o ^ with exponential 
decay at z — )• —00. 

The Evans function is now introduced as the following Wronskian: 

(4.3) D{X) = det (z+(0. A), . . . , Z+{0, A), Z'^^iO, A), . . . , ^"(0, A)) . 

Away from the origin X = 0, D vanishes at A with Re A > if and only if A is an eigenvalue 
of C, corresponding to existence of a solution Z{z, A) of CZ = XZ, decaying to at both 
z — > ±00. Indeed, the multiplicity of the root is equal to the multiplicity of the eigenvalue 
|GJH IGJ2| IMaZSl IZ4j . The meaning of the multiplicity of the root of D at embedded 
eigenvalue A = is less obvious, but is always greater than or equal to the order of the 
embedded eigenvalue |MaZ3| IZ4j . 

In agreement with |MaZ31 IZ4j , we define the Evans stability condition: 

(D) D has no root in {Re A > 0} except for A = 0, which is the root of multiplicity i. 

Note that under assumption (H3), the condition (D) is equivalent to D having precisely 
i zeros in {Re A > 0}. 

4.2 Type of the shock 

Define: 

£ = dimension of the unstable subspace of DG{V-^) 

+ dimension of the stable subspace of DG{V-^) — dimV, 

where diml/ = 6 is the dimension of the whole space. Then, the hyperbolic shock {V-, V-^) 
is defined to be: 

(i) of Lax type if £ = 1, 

(ii) of overcompressive type if £ > 1, 

(iii) of undercompressive type if £ < 1. 

li i = i > 1 OT i < i = 1, with i as in (H3), then the viscous shock V is defined to 
be of pure Lax, overcompressive, or undercompressive type, according to the hyperbolic 
classification just above. Otherwise, V is defined as mixed under- overcompressive type, 
|LZul IZHl IMaZ3t IZ4| . All the shocks considered in this paper appear to be of pure type. 
Indeed, though artificial examples are easily constructed |LZul [ZHj . we do not know of any 
physical example of a mixed-type shock. 



4.3 Linear and nonlinear stabihty 

Consider a planar viscoelastic shock for which the endstates V± lie in the region of strict 
hyperbolicity and profile {1^(-)} lies in the region for which i3o,j is uniformly elliptic. 
We have the following basic results relating the Evans condition (D) to stability. 
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Proposition 4.1 ( [MaZ3j ). A ssume (HO), (H2) and (H3). The Evans condition (D) is 
necessary and sufficient for the linearized stability Ci ^ of V, for all 1 < p < oo: 

||e*^/||L. <C(||/||ix + ||/|U.). 

Proposition 4.2 f |MaZ4l iRZj ). A ssume (HO), (H2), (H3) and (D). Then we have: 
(i) Stability. For any initial data V{-,0) with: 

Eo:= ||(l + |z|2)3/4(y(,o)-F)||^s«l 



sufficiently small, a solution V of (3.11) exists for all t > and: 

(4.4) 11(1 + \z\Y/\Vi;t) - V{- - stmnr. < CEo. 
(ii) Phctsc-ctsymptotic orbital stability. There exist Q^(t) (Xfid o^qq 

such that: 

(4.5) \\V{-,t) -V''^'^- - st)\\LP <CEo{l + t)-^^-^M/2 
and: 

(4.6) \a{t)-a^\<CEo{l + t)-^/^, \ait)\ < CEo{l + 1)-\ 
for all I < p < oo. 

Lemma 4.3. Assume (HO), (H2) and (D). If £ = i or i = 1, with i as in (D), then 
(H3) holds with the same value i. In particular, these conditions together imply nonlinear 
time- asymptotic orbital stability. 

Proof. The claim follows by the existence theory of |MaZ3j , relating the dimensions of stable 
and unstable manifolds of the rest points V± in the traveling- wave ODE, to the hyperbolic 
index i. Further |GZt IZHl IMaZ3j , stability condition (D) implies "maximal transversality" 
consistent with existence of a profile of the traveling- wave connection as a solution of the 
traveling-wave ODE (i.e. actual transversality), yielding (H3) with £ = i in the case i > 1, 
and (H3) with £ = 1. U 



Combining Proposition |4.2| with Lemma 4.3 we obtain: 



Theorem 4.4. For each of the planar systems considered in this paper, every viscous Lax, 
overcompressive, or undercompressive shock satisfying: 

(i) condition (H2) (noncharacteristicity), 

(ii) with endstates lying in the region of strict hyperbolicity of DG, 

(iii) with profile lying in the region of uniform ellipticity of Bq i, 

(iv) satisfying (D), 

is linearly and nonlinearly orbitally stable. 
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In particular, Propositions 4.1 4.2 



and Theorem 4.4 apply to profiles with a G such 



that the corresponding F of the form (3.1 ) is contained in a sufficiently small neighborhood 
of SO{2>). In the incompressible shear case, they apply to any profile with endstates a± ^ 0. 

The condition (H2) corresponds to noncharacteristicity of the shock, which holds gener- 
ically. It guarantees also exponential decay of the shock to its endstates |MaZ3l IZ4j , which 
is needed for efficient numerical approximation of the profile. 

Finally, we remark that strict hyperbolicity at V± is not necessary for existence of pro- 
files, but only to apply the basic stability framework developed in this section. When 
hyperbolicity fails, the corresponding endstate is unstable as a constant solution; however, 
this instability can be stabilized by convective effects if unstable modes are convected suffi- 
ciently rapidly into the shock zone; see Appendix [Cj This situation cannot occur for shear 
flows, for which all states are hyperbolic, but would be interesting to investigate in the 
compressible case. 



4.4 The integrated Evans condition 



Making the substitution V{z) = J^^V{y) dy and integrating the equations in (4.2) from 
— oo to z, we obtain after dropping the tilde notation: 

(4.7) XV = CV := BV" - QV' . 



We conclude for any A 7^ 0, satisfaction of (4.2) for a solution V decaying expon enti ally up 



to one derivative, implies that V{z) is also exponentially decaying and satisfies (4.7). 

Associated with C is an integrated Evans function D{\), which like D is analytically 
defined on the nonstable half-plane {Re A > 0}, through the construction sketched in 



section 4.1 In the Lax and overcompressive cases, the change to integrated coordinates 
has the effect of removing the zeros of D at the origin, making the Evans function easier to 
compute numerically and hence the stability condition easier to verify. 

Proposition 4.5 ( |ZHt IMaZ3| ). Assume (HO), (H2). Then the Evans condition (D) is 
equivalent to the following integrated Evans condition.- 

(i) for the Lax and overcompressive shock types: 

(D) the integrated Evans function D is nonvanishing on {Re A > 0}, 



(ii) for the undercompressive shock type: 

(D') the function D has on {Re A > 0} a single zero of multiplicity l + \i\, at X = 0. 

Note that the inclusion of term \l\ repairs an omission in |HLZj . for which £ = in the 
undercompressive case. Propositions |4.5[ |4.4[ [4^ and |4. 1| give together a simple and readily 
numerically evaluated test for stability of large- amplitude and or non-Lax- type waves. 
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4.5 Small-amplitude stability 



The following proposition gives a first nonlinear stability result for planar viscoelastic shocks, 
answering a conjecture posed in [AM ] for the shear wave case. 

Proposition 4.6 ( |HuZj ) . Assume (HO). Let Vq be a point of strict hyperbolicity of DG 
and let Aq be one of its eigenvalues, associated with a genuinely nonlinear characteristic 
field. Then there exists e > sufficiently small such that for any viscous shock V with speed 
s satisfying: 

\\y — Vb||L°° < e o,nd \s — Ao| < e, 

we have: 

(i) the shock is of Lax type, 

(ii) the Evans condition (D) holds, hence V is linearly and nonlinearly phase-asympto- 
tically orbitally stable. 



5 Existence of viscous shock profiles 

Let us now seek traveling waves connecting given endstates: 

y_ = y(_oo) = (a, 0), V+ = y(+oo) = (a+, b+). 



Indeed, by invariance of (2.1 ) under change in coordinate frame ^ i— t- ^-\-bot, we may without 
loss of generality assume that 6(— oo) = 0. 

Hereafter we restrict to the simpler (and apparently more physical) case of viscosity 
tensor Z2. The case Zi may be treated similarly. We note that the type and location of 
equilibria of the traveling wave ODE under are assumptions are independent of the choice 
of -Z, by the general results of |MaZ3] : see |BLZj for further discussion in the somewhat 
similar context of MHD. 



Writing the profile equation for (3.11) with (2.7) and ( |2.9| ), we obtain: 

■ib[,b'2,2bi,^ 



(5.1) -sa'-b' = 0, -sb'-DWoia)' 



03 



Note that s ^ for profiles satisfying the nonlinear stability conditions (namely, the end- 
states belonging to the strict hyperbolicity region of DG). For otherwise b' = and 
DWo{ay = hence M{a)a' = along the profile, contradicting the invertibility of M 
in the neighborhood of a(— oo). 

Now, substituting the first equation into the second, making the change of variable 
z I— >• sz, and defining cr = s^, we get the following reduced profile equation: 

(5.2) -aa' + DWo{a)' = { ^"'^^ 

V as 
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recognized as associated with the strictly parabohc gradient flux system in a alone: 

(01,02,203), 



(5.3) 
Note that rj{a) 



at + DWQ{a), 



as 



= ^ is the convex entropy for (5.2) as: 

V7?(o) • DlW{a) = Vq(a), q{a) = a ■ DaW{a) - W{a) 



Evidently, (5.2) may be written as a generalized gradient flow: 

|2 



(5.4) 



(o;,a'2,2o'3; 



a 



a 



03 



(o) = Wo(o) — {DWo{a) — aa) ■ o. 



where a = o(— 00). Making the change of variable z 1— )• z{z) where z solves the ODE: 
z'{z) = l/o3(z(z)), the system (5.4) becomes: 

(5.5) (01,02,203) = Va4>{a). 
We see that the function z 1— )• (j){a{z)) is non-decreasing: 

(5.6) ((/. o o)' = (V(A)o' = 03 diag (^1, 1, 



V0(o) (g) V(/) > 



in the admissible region 03 > 0. This is a simple instance of a more general fact concerning 
parabolic conservation laws possessing a viscosity-compatible strictly convex entropy [G] 
ICSll [U52| IBLZj . Moreover, the type of the shock connection of the original viscoelasticity 
equations is the same as the type for the reduced equations (5.3), which is in turn determined 
by the relative Morse index of the endstates/equilibria considered as critical points o: 

DWo{a) -aa- {DWo{a) - aa) = 

of (p. See also the general results and discussion of |MaZ3^ iBLZj . 
Finally, a straightforward calculation shows that: 



(5.7) 



S(/.(a) =sr]{V) - {q{V) + () + Vq{V) ■ (G(y) - G{V-) - s{V - y„)) 



for V = (o, b) with b = — s(o — a). 



where the inviscid flux G, entropy r] and entropy flux q are as in (3.3) and (3.4). The relation 
b = — s(o — a) is valid along the profile, and it follows by integrating the first equation in 
(5.1) from —00 to z. The vector C, = sDW{a)a — js^lap, which is independent of V, can 
be seen as an adjustment of the entropy flux q, naturally defined up to a constant. 
The quantity in the right hand side of (5.7) is related to the dissipative quantity: 

^^V) = -sr]{V) + q{V) 
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which decreases across any viscous profile connection lying within the region of strict hy- 
perbolicity of the reference hyperbolic-parabolic system and the region of strict convexity 
of its entropy rj (see |BLZ] and references therein): 

i;{V+)-i;{V-)<0. 



Indeed, by (5.7) and the Rankine-Hugoniot relations, it follows that: 
tl;{V-) = -s(p{a) and i;{V+) = -s0(a+). 



Thus, in view of ( |5.6| ), we conclude that in the present setting -0 is decreases across any 
viscous profile with positive speed s > 0, even one passing the elliptic region. This clarifies 
somewhat the role of (j) in the original system. 

5.1 The 3D compressible system 



Recalling ( |3.5| ), (5.5) becomes: 

a'l = (|ap — (T)ai — (|ap — a)ai, 
(5.8) 02 = (|ap — <T)a2 — (|a|^ — c)a2, 

2a'r^ = (|ap — 1 — cr)a3 — (|ap — 1 — a)as. 



As \a\ — )• oo, <j){a) ~ hence the phase portrait of (5.4) always possesses a minimum 
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or repellor. More, V0(a) ~ |apa points in the outward radial direction, and hence the 
index of this vector field on a suitably large ball is +1, and it must be equal to the sum of 



the indices of the equilibria (generically five - see Section 5.4 and Figure [2]), defined as the 



signs of the associated Jacobians sgndet(-D^VFo — fxld). The same argument shows that a 
sufficiently large ball is absorbing in backwards z, so that we can conclude that any orbit 
lying in the stable manifold of an equilibrium must connect in backward z to some other 
equilibrium possessing an unstable manifold. 

Further, when 03 = we have d^cj) = — — 1 — 5^)03, which is independent of 
(ai, 02). Since Vcj) ~ |apa as \a\ — )• +c«, it follows that the index of Vcp on a large half-ball: 
Br{0) n {03 > 0} equals +1 for (|a|^ — 1 — a)a3 > 0, and it equals for (|ap — 1 — 17)03 < 0. 
In the former case, the region Bji(0) D {03 > 0} is invariant in backward z and so we 
may conclude that any orbit lying in the stable manifold of an equilibrium in {03 > 0} 
must connect in backward z to some other equilibrium in {03 > 0} possessing an unstable 
manifold. 

5.2 The 2D incompressible shear case 



The incompressible case can be analyzed similarly as above, with (5.3) becoming: 

at + DWo{a)z = a^z, 
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where a 
model: 

(5.9) 



(ai,a2) and Wo{a) 



+ 5|op. That gives a 2 X 2 rotationahy symmetric 



a* + ((|ap + l)a). 



of a form at + {h{\a\)a)z = azz that has been much studied as a prototypical example of a 
system with rotational degeneracy [F]. Further, the counterpart of (5.4) reads: 

Ul2 



Va(/)(a) 



0(a) = Wo(a) 



a- 



{DWo{a) — aa) ■ a. 



Expanded in coordinate form, the profile ODE reads: 

|2 , 1 _\„ n„\2 



I 

0.2 



m- + 1 

|2 



a 



+ 1 



o-)a2 



ar + 1 

2 



+ 1 



cr)a2- 



Again, we find by an asymptotic development of (j) that the vector index of this ODE on a 
suitably large ball is +1, and there exists always at least one repellor, with index +1. In 
the generic case (see below), there are three nondegenerate equilibria, each of index ±1: 
one is of index +1 and of one index —1. 

Writing a = a(— oo) and a+ = a(+oo), the Rankine-Hugoniot relations for (5.9) are: 



(5.10) 



(|a+|^ + 1 - cr)a+ - (lap + 1 - a)a = 0. 



By rotation invariance, we may restrict our attention to a = (ai,0). We shall distinguish 
two cases. 



Case (i) a\ = 0. We find that there is a circle of solutions a+ to (5.10), given by 
|a+p = fj — 1, surrounding the rest state a+ = at the center. Along each radius of the 
circle, there is a viscous shock connection a(t, z) = p{z — at)e^^ solution to (5.9) whose norm 
p = |a| satisfies: 

(5.11) pt + {p^ + p). 



and connects p+ = p(+oo) = ^y a — 1 to p- 
associated parabolic equation to the flow: 



Pzz 

p{-oo) 



0. Note that (5.11) is also the 



(a? + l 



a ai 



aja. 



which is the counterpart of ( |5.4| ) for the Id incompressible model (3.19). When ai 
ai(— oo) = then ai+ = ai(+oo) = \/a — 1 and thus we obtain: 



/ 



{al 



a\)ai, 



which has the explicit solution: 
(5.12) ai{z) 



ai exp(— a^z) 
\/k + exp(— 2a|; 



with A; > 0. Note that this solution connects a(— oo) = \/ a — 1 to a(+oo) = 0; that is, the 
connection goes in opposite direction from the one sought. Setting now a(t, z) = ai {z—at)e^^ 
(with constant rotation angle 6) gives the traveling viscous shock solution to (|5.9|). 
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Remark 5.1. Though noncharacteristic when considered as one- dimensional solutions, as 
reflected by uniform exponential convergence to their endstates (see discussion below The- 
orem\4-4), such shocks are always characteristic with respect to the transverse (rotational) 



modes, which have characteristic speeds iby^|aP + 1 equal to ±\/cj = is 



Case (ii) a\ / 0. When 02+ = then (5.10) reduces to (af_|_ + l — (T)ai+ = {ai + l—a)ai 



corresponding to the associated scalar equation (5.11) 



When 02+ 7^ then the second equation in (5.10) becomes a\, + a 



whence, from the first equation: aj = \a\ = a — 1. That is, solutions with 02+ 7^ exist 
only if a is equal to the linearly degenerate characteristic speed, which has no profile. 

Thus we may without loss of generality restrict to the (at most) triples of possible rest 
states (a^*\o) with (to fix the ideas): 



a 



(3) 



J2) 



< < a\'' < 



(1) 



and ((aP)2 + l — cr)a\ < sufficiently small. Considering the equation (|5.1l[), we find that 



„(0 



the outermost rest points {a\ ,0) and {a\''',0) are connected to the innermost (a^^'',0) by 
a scalar (Id) shock profile. Indeed: 



„(i) 



(2) 



(aW)2 + l 

J2)x2 



(T < < 3(a 

,.(2)n2 



1 — o", so {a^i \ 0) is a saddle. 



(aY')'^ + 1 - o- < 3{aY')'^ + 1 - C7 < 0, so (of \ 0) is an attractor. 



3(af ^)2 + 1 - cj > {a^')'^ + 1 - o- > 0, so {aY',0) is a repellor. 



,,(3)^2 



(3) 



The phase portrait thus consists of a family of overcompressive profiles connecting {a\ ,0) 
and (a[^\o), bounded by Lax shocks between (a^^O) and {a^^\o), and between {a^^\o) 

and {a\ ,0), similarly as for the closely related "cubic model" studied, e.g., in [Fl|Br]. See 
Figure [1] for typical phase portraits computed numerically using MATLAB. 



5.3 The 2D compressible case 



Here (5.3) simplifies to: 



(5.13) 



a2,t + (|apa2) 



as,* + ((|a|^ - 1)03)2 



a2,z 
as 
as 



as 



where a = (02,03) and Rankine-Hugoniot relations are: 
(5.14) 



o+P — o')a2+ — (|ap — cr)a2 = 



(|a+| - 1 - a)as+ - {\a\ -l-a)as = 0, 
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(ay 



Figure 1: Typical phase portrait for the 2D shear case. In (a) we have ai = 1 and a = 11/4, 
and in (b) we have ai = 5 and a = 121/4. 



where a = (02,03) = a(— 00), 0+ = (a2+,a3+) = a(-|-oo) and 03+, 03 > 0. We shall 
distinguish two cases. 



Case (i) 02 = 0- Strict hyperbolicity of (3.24) enforces that 03 > l/\/3 and as / l/\/2. 



When a2+ = then (5.14) implies that: 

03+ = -^"3 ± ^Y^4(l + f7) - 3a|, 

with at most one physically feasible solution 03+ > 0. For every 03 the range of a, for 
which 03+ > and 03+ 7^ l/\/2 is: 

-,+00) \ jag + -^03 - -}. 



An associated Id traveling wave of the type (0, 03 (z)) must satisfy: 



' I / 3 
-0-03 + (03 



as) 



as 



If 02+ / 0, then the first equation in (5.14) implies |a+p = a, while by the second 
equation: 03+ = (1 + cj — a^ag, hence: 



a2+ = ±y a - [l + a - a^fal- 

We see that there is a pattern of at most four physically feasible equilibria, corresponding 
to two Id solutions plus two more symmetrically disposed about the 03 axis. There are at 
most five equilibria in total, counting a fifth possible infeasible radial solution with 03 < 0. 
Here, we are ignoring the line of nonphysical equilibria 03 = induced by the form of the 
viscosity tensor. See Figures [2] (a) (c) for a typical phase portrait computed numerically 
using MATLAB. 
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Figure 2: Typical three and five-equilibrium phase portraits for the 2D compressible case. 
The dark dashed lines bound the physically relevant region 03 > and the light dashed 



lines surround the region 7712 < where the system (3.17) loses hyperbolicity, see Section 
3.3.3 In (a) we have 0.2 = 0, 03 = 0.6, and a = 0.64, in (b) 02 = 0.2, 03 = 0.1, and o" = 4, 
in (c) a2 = 0, 03 = 2, cr = 4.84, and in (d) 02 = 0.2, 03 = 2, a = 5.29. 
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Case (ii) 02 / 0. Setting x = — |ap we see that (5.14) is solved by: 



a 



a 



a2A 



-03. 



.. . 1-^) 
(x+|ap— o") ' (x + |aP — 1 — cr 

Substituting into the definition of x and rearranging we obtain: 

{x + |q|^)(x + lap - af{x + jap - l-af =(|a|^ - cr)^(x + |a|^ - 1 - cr)^ai 
^^■^^^ +(x + |a|2-a)2(|a|2-l-cT)2ai, 



yielding the quintic: 

(5.16) y{y-af{y-l-af 



{\a\'-aY{y-l-aYai + {y-aY{\a 



l-ufal 



where y = |a+| . From the roots of (5.16) may be recovered through: 

|2 



(5.17) 



|a|" — a \a\" — 1 — a 
0.2+ = "2, ^3+ = — ; —03- 



y-a y-l-a 
Evidently, these solutions are not Id, as a2+/a3+ / ^2/03, unless y = |ap, in which case 



other roots of (5.16) are: a, 1 + a, contradicting (5.17) 



Recall that the nonphysical solutions with 03+ < are discarded and that the further 
condition of hyperbolicity of endstates is not necessary for existence of profiles, but is 
needed to apply the basic stability framework of Section |4j We shall discuss profiles with 
nonhyperbolic endstates in Appendix [Cj See Figures [2] (b) (d) for a typical phase portrait 
computed numerically using MATLAB. 

5.4 The 3D compressible case - continued 



In the full 3D compressible case, (5.3) can be written as: 



dt + (|a|^a) 



a3,t + ((lap - 1)03)2 



03, 

Q3,2 

03 



with a = (a, 03) and a = (oi, 02). The corresponding Rankine-Hugoniot relations read: 

(|a+p — a)d^ — (jap — a)a = 
(|a+p - 1 - cr)a3+ - (|ap - 1 - £7)03 = 



(5.18) 



where a = (0,03) = a(— 00) and a+ = (a+,a3+) = a(+oo) and 03^,03 > 0. Also, by 
invariance with respect to rotations in the a plane, we will assume that ai = 0, without 
loss of generality. 
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Case (i) 02 = 0. Here, the phase portrait can be easily deduced from that in (i) Section 
5.3 of the 2d compressible case. That is, there is at most one physically feasible Id profile 



connecting to rest point a+ = (0,0, (—as + \/ 4(1 + a) — 3a|)/2), and a ring of rest points 



a+ = (r cos 0,r sin 0,03+) with 03+ = (1 + a — 03)03 and r = cr — a|_|_. Again, we are 
ignoring possible equilibria in the nonphysical plane 03 = 0. 



Case (ii) 02 7^ 0. This case includes the 2d portrait of case (ii) in Section 5.3 for the 
2D compressible case when oi = 0. 

If ai+ 7^ 0, then |a+P = o", and so |a2p = o" (in view of the second equation in (5.18)). 
Thus, except in this degenerate case, the set of equilibria is only that of the planar case 
already treated. The types of shock connections may be different than in Section 5.3, and 
profiles may go out of plane to yield new connections. 

The above situation is quite reminiscent of the case of MHD [BLZj . In particular, if 02 
is varied slightly from the rotationally symmetric situation 02 = 0, then one may conclude 
by persistence of invariant sets as in |FSj that "Alfven"-type profiles must arise in the 
rotationally degenerate characteristic field. 



6 Numerical stability analysis 

In this section, we describe the numerical Evans function method, based on the numerical 
approximation using the polar-coordinate algorithm developed in |HuZj : see also |BHRZ| 
IHLZ| |HLyZ[ IBHZ] . Since the Evans function is analytic in the region {ReX > 0} of interest, 
we can numerically compute its winding number around a large semicircle .6(0, R)r]{ReX > 
0}, enclosing all possible nonstable roots. This allows us to determine stability through the 
Evans condition (D); alternatively, as we shall do here, through its integrated version (D) 
(resp., (£)')). In the case of instability, one may go further to locate the roots and study 
stability and bifurcation boundaries as model parameters are varied. This approach was 
introduced in basic form by Evans and Feroe |EFj and it has since been elaborated and 
greatly generalized. For applications to successively more complicated systems, see for 
example [PSWl RSl IB^ lB?Zl iBDGl iHuZl iHLZl |HLyZ[ iBHZllBLZ] . 



6.1 The Evans systems 



Linearizing about a traveling wave solution (a, b) = (oi, 02, 03, 61, 62, 63) of (3.14), we obtain 
the eigenvalue problem: 

Xttj — sa'j — b'j = for j = 1, 2 
- sa's - 63 = 
Xbj — sb'j — (lapOj + 2(a • a)aj)' = {b'j/as — asbj/a,'^)' 
A63 - S63 - ((|a|^ - 1)03 + 2(a • 0)03)' = 2(63/03 - ashyal)'. 
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We make the substitution di{z) = J^^ai{y) dy and bi{z) = J^^bi{y) dy into (6.1) and 
then integrate from — oo to z to obtain, after dropping the tilde notation: 



Xttj — sa'j — b'j = 



(6.2) 



Aas — SO3 — &3 = 
Xbj — sb'j — {\a\'^a'j + 2{a ■ a')aj) = bj /a^ — a'^b'j/a'^ 
A63 - S63 - ((lap - 1)03 + 2(a • a')a3) = 2{bl/a3 - 463/0!). 



6.1.1 The 3D compressible case 



In the full 3D case (6.2) may be written as a first order system Z' = A{z, X)Z, with 

Z = (61,01,4,62,02,02,63,03,03)^ 

and 

(6.3) A{z,X): 
where: 

OA -s 



(6.4) 



B{z,X) C{z,X) 
D{z,X) E{z,X) 



Biz,X) 







-Aa3 , A + a3(|ap-s2 + 24) 
A03 



C{z,X) 



s s 
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D{z,X) 
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(6.5) 






s 





A 







—s 



X + asilal^ - + 2al) 



a2al 



A 







-b'2 + 20203 

503 



1 



-Aa3 Aa3 a^(|ap - 1 - ^2 + 2a^) - 26^3 + 2a3A 
2s ^ 2sS3 



6.1.2 The 2D incompressible shear case 



For (3.15), the same procedure as in (6.1) yields: 

Xa — sa — 6' = 0, 



(6.6) 



Xb - sb' - ((1 + |a|2)a + 2(a • a)a)' = b" . 



Substituting ai{z) = f^^ai{x) dx, bi{z) = f^^bi{y) dy, into (6.6) and integrating from 
— c« to z we obtain, after dropping the tilde notation: 



(6.7) 



Aoj — sa[ — 6'j = 0, for i = 1, 2 
Xbi - sb'i - ((1 + |a|2)a- + 2(a • a')ai) = 6'/. 

Let Z = (ai, 61, 6']^, 02, 62, 62)"^. Then (6.7) may be written as (4.2), where: 
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Figure 3: (a) Traveling wave profile V for the shear case with parameters values ai = 1, 
02 = 0, and s = 1.8547 corresponding to a Lax shock connecting endstates (1,0) an (0.8, 0). 
(b) The image of the semicircle under Evans function D. 



(a) 



(6) -'-^ 



Figure 4: (a) Traveling wave profile V for the shear case with parameter values ai = 1, 
02 = 0, s = 1.8547 corresponding to an overcompressive wave connecting endstates (0.8,0) 
and ( — 1.8,0). (b) The image of the semicircle under Evans function D. 



6.1.3 The 2D compressible case 



With j = 2 in (6.2) and using b'j = Xaj — sa'j, (6.2) can be equivalently written as (4.2) 
with Z = (62, «2, ^3) fl's)"^ and A{z,X) = E{z,X) given in (6.5) where a = (02,03). 
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6.1.4 Transverse equations 



Consider now a 2D compressible solution as a solution of the full 3D system (3.14). We 



find that the integrated eigenvalue equations (6.2) decouple into the 2D equations plus the 



transverse system, obtained from the equations corresponding to j 
putting ai = bi = 0: 



1 in system (6.2) after 



Aai 



sa'i 



b'l = 0, 



Xbi — sb\ — |apa' 



as 



The system (6.8) has the form of as (4.2) with Z = (5i, ai, a'l) and A{z^ A) = A) given 
in (6.4), after putting oi = 0. 



6.2 Approximation of the profile and of the Evans function 



Following |BHRZ1 IHLZ| . we approximate the traveling wave profile using one of MAT- 
LAB's boundary- value solvers bvp4c |SGT| . bvp5c |KLj . or bvp6c |HM] . These are adap- 
tive Lobatto quadrature schemes that can be interchanged for our purposes; for rigorous 
error /convergence bounds for such algorithms, see e.g. |BeH IBe2j . The calculations are 
performed on a finite computational domain [— L,L], where the values of approximate plus 
and minus spatial infinity L are determined experimentally by the requirement that the 
absolute error |y(ibL) — V±| < TOL be within a prescribed tolerance, say TOL = 10~^. 
Here V and V± are the profile and limiting endstates as defined in Section 



4.1 



Using now the notation of section l4?Tl define for z > and z < 0, respectively: 



Z+(z,A) = Z+(z,A) A...AZ+(z,A) and {z,\) = Z^^^{z,\) h . . . h Z^{z,\), 
so that the Evans function is given by: 

(6.9) L»(A) = Z+(0,A) A2"(0,A). 

Since for L > large, P±(±L, A) approximately equals Id, we obtain: 

(6.10) Zf{±L, A) ~ e^-^i^^^^Zf (A), 

where the "+" sign is taken with indices i = 1 . . .k and the "— " sign with i = k + 
1 . . . N. Recall that 5(A) = span{Z^ {X)}i=i jt is the stable space of ^+(A), while 1{{X) = 
span{Zlf^ {X)}i=k+i..N represents the unstable space of ^_(A). 

The analytic bases {Zf{X)}i are obtained by the following procedure |HuZl [Z2]. First, 
one computes the eigenprojections V+{X) and 'P-(A) onto, respectively, the space 5(A) and 
U{X). This can be done by setting: 



V± = R±{L^R±) 'l±. 
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where R±{X) and L±(X) are matrices consisting of any orthonormal right and left bases of 
5(A) (when with subscript "+") and (when with subscript "-"). Then, a standard 

result in matrix perturbation theory [K] states that the analytic in A bases ^^(A) can be 
prescribed constructively as the solution of Kato's ODE: 



{Z^y = {V'^V± - V±V'±)Z^, Z±(Ao) = ii±(Ao), 

where ' denotes the differentiation with respect to A. This prescription is also minimal in 
the sense that V±R'^ = 0. 

We now continue the construction of the approximate Evans function. As a consequence 
of (6.10), for large L we set: 

Z+{L, A) ~ ^+,(L, A) := etr('^+Wl^(^))^Z+(A) A ... A Z+(A), 
Z-{-L, A) ~ 2-p,(-L, A) := e't^^"^- Wl"w)^Z-_^i(A) A ... A Z^(A). 

The objective is now to trace the evolution of the differential form Z'^pp[-, A) backward in 
z, and the evolution of Z^p(-,A) forward in z, starting from, respectively, the initial data 
Z^pp{L, A) and Z~pp{—L, A), and according to the system as in (4.2): 

(6.11) Z'Jz, A) = A{z, X)Zapp{z, A). 



The numerical approximation of D{X) in (6.9) is then recovered through: 
(6.12) D{\) ~ D,pp{\) := Z^ppiO, A) A ^"^^(0, A). 



To solve (6.11) for Z^pp we use the polar-coordinate method described in |HuZj . which 
encodes Z^pp as product of a complex scalar and the exterior product Q,^ of an orthonor- 
mal basis {ujf} of S or, respectively, an orthonormal basis {w^"} of ^: 



Zt„p{z,X)=r^{z,X)n^{z,X), 



uj{ A ... Aoj^, 



The above quantities evolve by some implementation (e.g. Drury's method below) of 
continuous orthogonalization, where the "radius" r satisfies a scalar ODE slaved to 0, 



related to Abel's formula for evolution of a full Wronskian. Namely (6.11 ), is equivalent to: 

Q'(z, A) = (idN - nn*'^A{z, X)n{z, x) 

r'{z, A) = tiin* A{z, X)n^ • r{z, A) 



(6.13) 



and we recover, in view of (6.12): 



DappiX) = r+(0, A)r-(0, A) • l^+(0. A) A ^-(0, A), 



see fHuZj ESI EH] for further details. The rationale for solving the system (6.13) for the de- 
composition of Zapp, rather than the original (6.11 ) is that the imposition of orthonormality 
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on 17 prevents the collapse of the various columns (solutions) onto a single fastest-growing 
mode, as would otherwise be the case. For a discussion of this and other numerical issues 
connected with the polar coordinate method, see |HuZ| [Z3] . 



The calculations of (6.13) for individual A are carried out using MATLAB's ode45 
routine, an adaptive 4th-order Runge-Kutta-Fehlberg method (RKF45) with excellent ac- 
curacy and automatic error control. Typical runs involved roughly 60 mesh points per side, 
with error tolerance set to AbsTol = le-8 and RelTol = le-6. To produce analytically 
varying Evans function output, the initializing bases {Z^} are chosen analytically using 
Kato's ODE [GZl iHuZl iBrZl IBHZ] . Numerical integration of Kato's ODE is carried out 
using a second-order algorithm introduced in \Z2\ IZ3] . 

6.3 Winding number computation 

Recall that the Evans condition amounts to checkin g for the existence of unstable zeros of 



the integrated Evans function D, described in section 4.4 We first observe (Proposition 6.10 
|HLZj ). that for shock profiles of the hyperbolic-parabolic systems of the type we consider, 
there holds: 

(6.14) lim — ^-2 = C uniformly on Re A > 0, 

|A|-s>oo gOVA 

with constants a and C ^ 0. When D is initialized in the standard way on the real axis. 



so that -D(A) = -D(A), a and C are necessarily real. The knowledge that limit in (6.14) 
exists allows actually to determine a and C by curve fitting of logZ)(A) = logC + oA^ 
with respect to A^/^, for large |A|. 

One further determines the radius R > so that: 

l)(A)/0 for |A| > i? and i?e A > 0, 

by taking i? to be a value for which the relative error between D{X) and Ce°^ becomes 
less than 0.2 on the entire semicircle: 



= d(^B{0,R)n{Re A > 0}j, 

indicating sufficient convergence to ensure nonvanishing. For many parameter combinations, 
R = 2 was sufficiently large. Alternatively, we could use energy estimates or direct tracking 
bounds as in |HLZj and |HLyZ| , respectively, to eliminate the possibility of eigenvalues of 
sufficiently high frequency. However, we have found the convergence study to be much more 
efficient in practice; see |HLyZ| . 

We now compute the winding number I{R) of the image curve D^Sj^) with respect to 
0, which equals the degree of the 2d vector field given by D in the interior region of Sr. 
Since the index of any nondegenerate zero of a holomorphic function is +1, the condition 

I{R) = 
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is hence equivalent to D having no zeros in the open interior of the curve Sr. Since all the 
shocks considered here are of Lax or overcompressive type, condition I{R) = is equivalent 
to the Evans stability condition {D). 

The winding number I{R) is now computed by varying values of A along 20 points of 
the contour 5, with mesh size taken quadratic in modulus to concentrate sample points 
near the origin where angles change more quickly, and summing the resulting changes in 
arg(Z)(A)), using 91og-D(A) = argl)(A)(mod27r). To ensure winding number accuracy, we 
test a posteriori that the change in D for each step is less than 0.2, and add mesh points as 
necessary to achieve this. (Recall, by Rouche's Theorem, that accuracy is preserved so long 
as relative variation of D along each mesh interval is < 1.0.) In Tables [i] and [l] we give the 
radius of the domain contour, the number of mesh points, the relative error for change in 
argument of D{\) between steps, and the numerical approximation of spatial infinity ztL. 



6.4 Results of numerical experiments 

In our numerical study, we sampled from a broad range of parameters and checked stability 
of the resulting Lax and over-compressive profiles whenever their endstates fell into the 
hyperbolic region as required by our stability framework. We did not find any undercom- 
pressive profiles for the model considered here, either in the incompressible shear or the 
compressible case, nor did Antman and Malek-Madani find undercompressive profiles in 
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Table 1: Table demonstrating contour radius, number of mesh points, relative error, and 
spatial domain for the incompressible case. 
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Table 2: Table demonstrating contour radius, number of mesh points, relative error, and 
spatial domain. The data on the left side corresponds to the compressible 2D system, and 
that on the right to the transverse system. 
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their investigations of the incompressible shear case |AMj . As shown in Section |Dj under- 
compressive connections cannot occur in the incompressible shear case for any choice of 
potential. However, we do not see why they could not occur for other choices of elastic 
potential in the compressible case. 

All our computations yielded zero winding number, consistent with stability. All to- 
gether, our study consisted of over 8,000 Evans function computations. The following 
parameter combinations were examined for Evans stability. 

The 2D incompressible shear case. The following parameter combinations yielded 
Evans function output with winding number zero, consistent with stability: 

(a, s) E {0.2 : 0.2 : 5} x {0.2 : 0.2 : 7}. 

The 2D compressible case. In the compressible 2D case and the transverse case following, 
we computed the Evans function for the stated parameter combinations whenever the profile 
end-states did not lie in the elliptic region. For 02 / 0, we restricted our attention to the 



profiles connecting rest-points corresponding to solutions of (5.16) in the interval [—50, 50]. 

We computed the Evans function for all 2 point configurations (Lax connections) coming 
from the following parameter combinations. All computations yielded zero winding number, 
consistent with stability: 

(02, as, s) e {0} X {1.1 : 0.5 : 25.6} x {0.5 : 0.5 : 20}. 

(02, as, s) G {0.1 : 0.4 : 25.3} x {0.2 : 0.4 : 25.4} x {0.3 : 0.4 : 25.5}. 

For the following parameter combinations, we investigated the 4 point configurations com- 
puting all Lax connections, and 5 overcompressive connections passing through evenly 
spaced points on the segment in phase space connecting the saddle points: 

(a2,as,s) G {0} x {0.03 : 0.07 : 1.03} x {0.05 : 0.07 : 1.05}, 

(a2, as, s) = {(0.08, 0.59, 0.75), (0.08, 0.87, 0.82), (0.22, 0.66, 0.75)}. 

The transverse case. We computed the transverse Evans function for the two point 
configurations (Lax connections) for the following parameter combinations: 

(a2, as, s) G {0.1 : 0.4 : 25.3} x {0.2 : 0.4 : 25.4} x {0.3 : 0.4 : 25.5}. 

In addition we examined the 4 point configuration corresponding to (a2, as, s) = (0.1, 0.8, 0.8) 
computing the Evans function for the 4 Lax connections and for 5 overcompressive connec- 
tions passing through points evenly spaced along the line in phase space between the two 
saddle points. 

6.4.1 Numerical performance 

The Evans function computations for the most part worked reliably and well, showing 
performance comparable to that seen in previous studies for gas dynamics |HLZ1 HLyZ| 
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and MHD |BHZ1 IBLZj . A typical winding number computation for a single profile took 
approximately 30 seconds and computation of the profile approximately 5 seconds. 

As expected, performance degraded catastrophically in various boundary situations: the 
small-amplitude limit as — )• 0; the characteristic limit as one or more characteristic 
speeds approach the shock speed; the large-amplitude limit as |a-i-| approach infinity or 
as approaches the physical (infinite compression) boundary as = 0; and the elliptic limit 
as one or both endstates a± approach the elliptic region where characteristic speeds are 
complex. For discussion of causes of and (partial) cures for these numerical issues, see, e.g., 
|HLZ1 IBHZt IBLZl IZ3| . In the present study, such boundary cases were omitted. 

7 Discussion and open problems 

In this paper, we have obtained the first analytical stability results for viscoelastic shock 
waves, stability of small-amplitude Lax shocks, and set up a theoretical framework for 
future numerical and analytical studies of shock waves of essentially arbitrary viscoelastic 



models. A large-scale numerical Evans study for the canonical model (2.7) yielded a result 
of numerical stability for each of the more than 8,000 profiles tested, of both classical Lax 
and nonclassical overcompressive type, and with amplitudes varying from near zero to 50. 

Interesting problems for the future are the treatment of more realistic potentials with 
physically correct asymptotic behavior, systematic numerical and asymptotic investigation 



across parameters as in [HLZ] HLyZ , IBHZl IBLZ| , and the treatment of phase transitional 



elasticity by incorporation of dispersive surface energy terms. 

A Appendix: General facts 

Though the investigations of this paper were carried out for special choices of VF, Z, the 
methods we use apply to much more general choices. With an eye toward future work, we 
collect in this appendix the information needed to carry out such extensions. 

A.l General elastic potential 



Theorem A.l. Let W : M.^^^ — M_|_ satisfy (2.4) and (2.6). Then there exists a scalar 
function a : — > M+, such that: 



W{F) = a{\F\^, \FF^ |^detF). 

The derivative DW{F) e R^''^, wherever defined at F e R^""^ (so that dAW{F) = DW{F) : 
A), is given by: 

DW{F) = Vo-(|Fp, |FF'^|2,detF) • (^2F, 4FF'^F, cof . 

IfW(ld) = and W is in a neighborhood of S0{2>), then: 

DW{ld) = 0, D'^W{ld) : A = A(trA)Id + /isymA G M3^^ 
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with the convention d\_^ y^M^(Id) = {D'^Wild) : A) : Ai and the Lame constants X and fi: 

A = VV(3, 3, 1) : ((2, 4, 1) ® (2, 4, 1)) , = Va(3, 3, 1) • (0, 8, -2) 
satisfying: fi > and 3X + fi > 0. 

Proof. According to the representation theorem jTNj . every frame invariant and isotropic 
W depends only on the principal invariants of the left Cauchy deformation tensor FF^ , that 
is W{F) = a{tT{FF^),ti cof (FF^), det(FF^)). Since tr cof Q = ^(tr - ^tr (Q^), the 
claim on the form of W follows directly. 

The formula for derivative DW{F) follows from: 

OaIFI^ = 2F:A, OaIFF^I'^ = AFF^F : A, a^detF = cof F : A, 

where for the last expression we used det(F + Q) = det F + F : cof Q + Q : cof F + det Q, 
valid for 3x3 matrices F, Q. 

The vanishing of DW{ld) is clear since W is minimized at Id. The formula for D^W{Id) 
follows by chain rule and Lemma A. 2 Further, notice that d'^WA,Ai^d) > 0, which reads: 

(A.l) A|tr + ^| sym > G R^''^. 



Evaluating (A.l) first at a traceless A and then at yl = Id we see that: 
(A. 2) /i>0, 3A + ;U>0. 



To prove that (A. 2 ) implies (A.l), write sym A a the sum of orthogonal matrices: sym A 
diag(oii, a22, 033) + B. Then: |sym = ^^^^ al + \B\^, so: 



33 3 

|2 



A|tr A\^ + fi\ sym A\'^ = A(^ au)^ + fi^al + fi\B\^ > (A + /u/3)(^ au)^ + fi\B\ 

1=1 i=l i=l 

which ends the proof. ■ 

For the behavior of W close to the energy well 50(3) it is important to know the 
derivatives of at i? G SO{3). It follows by frame invariance that: 

dlFu...RF„W{R) = a^,,...^„VF(Id) VFi . . . F„ G r3x3 ^ g^^^^y 
Hence, it suffices to find the derivatives of W at Id. Direct calculation yields the following: 

Lemma A.2. For F G R^''^, let a{F) = |F|2, /3(F) = jFF'^p and 7(F) = detF. Then, 
for any Ai,A G M^^^ we have: 

(a,/3,7)(Id) = (3,3,1) 
9Ai(a,/3,7)(Id)= (ld:^i)(2,4,l) 

dA,,Ai^, 7)(Id) = {a : ^1) (2, 4, 1) + (symA : Ai) (0, 8, -2) 

5ii,A,A(«' 7)(Id) = - (cofyl : ^1) (0, 8, -2) + ((cof^ + AA'^ + 2A^) : ^1) (0, 8, 0) 

4iAA,A(«,/3,7)(Id) = sUaA^'A) : ^1) (0,8,0) 
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For F as in ( |3.1| ), we have: 

a(F) = 2 + |a|2, /3(F) = 2 + |o|^ + 2(|o|2 - ai), 7(^) = «3. 

Hence, and without loss of generahty, the reduced function W{a) must be of the form 
W{a) = a(|ap,a3), for some : — ^ 1^. 

In the incompressible shear case, it reduces to: 

(A.3) W{a) = a{\a\^) = a{\a\'',l), 

leading to a profile equation agreeing with that of the 2x2 rotationally symmetric model: 
(A.4) at + i2Va{\a\'^)a), = a,,. 

This clarifies and puts in a more familiar context the investigations of Antman and Malek- 
Madani [AMj on existence of viscous profiles for the 2D incompressible shear model. 

A. 2 General viscous stress tensor 

We do not have a complete categorization of possible Z satisfying (i)-(iii) corresponding 
to that of Theorem |A.l for the elastic energy density W. However, we note the related 



discussion of Antman [A] for the class of systems of strain-rate type: 

ZiF,Q) = FSiC,D), 

where S* is a symmetric dissipation tensor depending on the the metric C = F'^F and on 
its time derivative D = F^Q + Q^F = 2sym.(F'^ Q) . These automatically satisfy (i) and 

(ii) because sym{{RF)'^ {RK F + RQ)) = sym(F'^Q + F^KF) = sym(F^Q). Thus only 

(iii) need be checked, in the form: 

(A.5) S{C,D):D>Q. 



Both of the examples (2.9) are of this type. For Zi = FSi, we have: 



Si{C,D) = sym{F^Q) = ^D, 



which evidently satisfies the strict version of inequality (A.5): 
(A.6) S{C, D):D> -f\D\^, 7 > 0. 

For Z2 = FS2, we have S2 = {detF)F-hym{QF'^)F-^'^ = l{det Cy/^C'^DC-^, hence: 

82(0, D) : D = ^{det Cy/^C-^CtC-^ : D = ^(det C7)i/2|C7-i/2^c'-i/2|2 > ^|^|2^ 
where 7 depends on C, but is uniform for F bounded and det F bounded away from zero. 



In [A], Antman proposes as a sufficient condition for (A.6), that S be monotone in D 
in the sense that : A > 7|Ap, 7 > 0, for A symmetric. This implies (A.6) under the 



additional assumption S'(-,0) = (i.e. viscous force vanishes at zero velocity), by: 

S{C,D) : D - S{C,0) : D = / —{C,eD)D:Dde>j\D\'^. 

Jo OD 

It is easily checked that monotonicity holds for both of the choices -Zi, Z2- 
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B Appendix: Phase-transitional elasticity 



Another interesting direction for future investigations is the phase-transitional case, as we 
now briefly discuss. A typical model, as described in [FPj . has the form: 



W{F) 



C4 



where: 



C± 



FlF± 



1 
1 ±e 
±e 1 + £2 



10 

1 ±e 
1 



Evidently, W is minimized among planar deformation gradients F at the two equilibria F± . 



A particularly interesting class of solutions to (2.1) are stationary phase-transitional 
shocks connecting the two equilibria {b±,a±), a± = (0,ibe, 1), that is, zero-speed shocks 
compatible with the Rankine-Hugoniot condition for ( |3.11 ) with (3.3): 6+ = DW{a+) = 
DW{a ). Similarly as in section [s] and as in Id case treated by Slemrod |S1], such connec- 
tions do not exist under the viscoelastic effects alone, but their existence requires also the 



inclusion of third-order surface energy terms as in Section 2.3 



Solving the system ( |5.1[ ) augmented by the term 7div £ coming from the surface energy 
£q as in (2.11), we obtain (for s = 0) that b = const and DW{a)z = "ya-zzz- Take the capil- 
larity coefficient 7 7^ and assume that DW{a±) = (the end-states in the potential well). 
Integrating, we obtain a harmonic oscillator equation: DW{a) = "jazz, with Hamiltonian: 



H{a,az) = Wia) 



-I 
2' 



const. 



Hence, a question is whether the connected component of the level set of H, containing 
(a_,0), contains also (a+,0). Contrary to the Id case in [SI], this is only a necessary 
condition for profile's existence. Existence or nonexistence of such connections would be an 
interesting question for further analytical and numerical investigation. 

A second question would be to determine stability of such stationary transitions, should 
they exist. We conjecture that, similarly as in |Z8j for the Id case, the spectral stability fol- 
lows automatically by energy-considerations. Stability of nonstationary phase-transitional 
profiles has not been treated even in the Id case, and would be another interesting problem. 

Finally, we point out that the equations with surface energy do not fit the stability theory 
of Section [4j since they are third- and not second-order. However, we expect that the basic 
methods should still apply, after suitable modifications. It would be very interesting, for 
the sake of this and other applications involving dispersive phenomena (for example, the 
Hall effect in MHD), to carry out a complete analysis extending the nonlinear stability 
framework to this higher-order case. 
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C Appendix: Nonhyperbolic endstates 



Consider the ID compressible equations ( |3.20[ ) with the associated profile equation, obtained 
by setting ai = 02 = in (5.8): 



(C.l) 



203 = (a| — a3 — eras) — (a^ 



a 



era), 



a := 03- 



As noticed in section 3.3.4, strict hyperbolicity of (3.20) corresponds to las] > l/VS. On 



the other hand, when a > 03+ > 0, + (t)/3 < a < ^/l + a and a = a'^_^_ + 03+0 + > 



0, the equation (C.l) posseses a solution as{z), decreasing from its unstable equilibrium 

03 (+00). This solution corresponds to a viscous shock 
03,22 and the 2d viscous system: 



a = a3(— 00) to the stable one 03+ 
profile of the associated scalar equation 03,* + {a^ — 03) 



(C.2) 



«3,t 



33,2 



0, 



53,t 



53,22- 



Thus, there exist shock profiles of (3.20) for which the left end-state is hyperbolic, with 



one characteristic greater than the shock speed s and the other smaller than s, but the 
right end-state is not hyperbolic as it has two pure imaginary characteristics with real parts 
c+ < s. This is a "complex Lax shock" of the type considered in tAMPZ| IOZ| 



Writing the system (C.2) in the operator form (03, bs)J = Q ■ (03, 63)^ where: 



Q 



(3ai - l)d. 



dz 



we determine the spectrum of the constant solution (a+, 6+) through the dispersion relation: 



= det 



Consequently: 



-(3a2 



m 



A 

3+ " 



l)ik 



—ik 
A + P 



A^ + k'^X + (3aL - l)k 



-k'^±k 



1 



for ~ 0, 



yielding a maximal growth rate e^*' with r = Re Amax(^) ~ (1 — 3a3_,_)/2. As discussed in 
[AMPZllOZj . for speed |s| sufficiently large compared to the growth rate r, the shock profile 
can be seen to be stable, despite instability of its right end-state as a constant solution. 
It also can be shown, by weighted coordinate techniques as in |Satt ILRTZj . that Evans 
stability together with such a convection vs. growth condition on the essential spectrum, 
implies linearized and nonlinear stability for perturbations that are exponentially localized 
on the half-line z > 0. Another direction for further investigation might be to understand 
whether there are interesting physical phenomena corresponding to shocks of this type. 

In contrast with the situations considered in [AMPZl lOZj , for which profiles associated 
with complex Lax shocks are oscillatory at the complex end, the profiles here are of ordinary 
monotone type. Other complex Lax connections, genuinely two-dimensional, may be seen 
in Fig. ^h). Provided that all characteristics are incoming on the complex side z — )• +00, 
nonlinear stability may again be established by weighted norm methods assuming Evans 
stability plus an appropriate convection vs. growth condition as described in |AMPZ| IOZ| . 



38 



D Appendix: Nonexistence of undercompressive profiles 



We show that undercompressive profiles do not occur for general shear models as in (|A.3) 



The same argument implies that undercompressive profiles cannot occur also in the 2D 
compressible case, for the special class of potentials W{a) depending only on |ap. 



Note first that every undercompressive profile {a{z), b{z)) of (3.15), with speed s, induces 
an undercompressive profile of (5.9) given by z i— )■ a{sz) with speed a = s^. This statement 
follows from a more general observation in |MaZ3l IBLZ| relating the type of inviscid shocks 



to the connection number for the traveling wave in the reduced ODE (5.2). Alternatively 



the same can be checked directly: if the sum of the number of eigenvalues A for (3.15) at 



a(— oo) with A > s, and the number of eigenvalues fi at a(+oo) with < s, is less than 5 
(the dimension of the system increased by 1), then the sum of the number of eigenvalues of 
(5.9) A^ > o" and the number of eigenvalues /x^ < o" is less than 3, when s > 0. When s < 



we likewise have that the number of < o" plus the number of /i > o", is less than 3. 



We shall prove that the system: 



(D.l) 



at + {h{\a\^)a) 



generalizing (5.9) and (A.4), where h = 2Vcj, admits no undercompressive shocks. 



After integrating in z, the traveling wave ODE for ( |D.l ) reads: 
(D.2) a' = -aa + h{\a\'^)a - ( - cra_ + /i(|a_p)a_), 

where a is the speed of the shock and a_ = (ai_,a2-) = a(— oo) is its left end-state. Note 
that all possible right end-states a+ = (ai+,a2+) which can be connected to a_, satisfying 
hence the Rankine-Hugoniot condition: 



(D.3) 



/i(|a- 



aa- 



/l(|a^ 



aa+, 



must lie on the same line through the origin, due to rotational invariance of the system (D.l ) 
We may, without loss of generality, assume it to be the ai axis, so that: a2- = 02+ = 0. 



The gradient of the right hand side in (D.2): 



D{h{\ 



a a 



aa 



) = {H\c 



cr)ld + 2V/i(|a| 



has two eigenvalues: Ai(a) = h{\a\'^)—a+2Vh{\a\'^)\a\'^, with the radial direction eigenvector 
a, and A2(o) = /i(|ap) — a with the eigenvector a-*- in the transverse (rotational) direction. 

Observe further that an undercompressive shock profile must necessarily be a saddle to 
saddle connection, and that any saddle point has one of its invariant manifolds (the one 
corresponding to Ai) confined to the oi axis. We find that the profile must either lie entirely 
on the ai axis, or else entirely off the axis, leaving o_ and entering a+ along the transverse 
direction (orthogonal to the axis). We now distinguish 3 cases: 

Case (i) ai+ai_ 

= h{al_ 



A2(a_ 



> 0. In this situation, in view of (|D.3|), the transverse eigenvalues 
a and A2(a+) 



1+ 



) — a must also have a common sign. Thus 
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the profiles may only leave or only enter along the transverse directions, contradicting the 
assumed behavior. 

Case (ii) ai+ai < 0, and the profile connection is radial. Without loss of 
generality, assume that ai_ > and ai-(_ > so that a'^{z) < along the whole profile. In 



particular, ai{zo) < at zq where ai{zQ) = 0. By (D.2) it follows that: 



/i(a^_)ai_ — (Tai_ = —a'i{z()) > 0, 

hence h{a1_)—a > 0. This means that the transverse eigenvalue corresponds to the unstable 
direction, contradicting the profile being radial and a_ being a saddle equilibrium point. 

Case (iii) ai+ai < 0, and the profile connection is transverse. In this situation 
A2(a-) > and A2(a+) < 0. Observe that on the circle |ap = af_, we have: 

a' = (/i(a^_) — o") (a — a_) = A2(a_)(a — 



and thus the exterior of this circle is a positively invariant set. Likewise, by (D.2) and (D.3), 
on the circle |ap = af_^_ there holds: 

a = (/i(ai+) - cr) (a - a+) = A2(a+)(a - a+), 

and hence the exterior of this circle is a negatively invariant set. Consequently, at that of 
the two saddle points a_ and a+ which has larger norm, the invariant curve tangent to the 
transverse direction must lie, for all times, outside the corresponding circle. This is clearly 
a contradiction and ends the proof, as the case ai+ = — ai_ is ruled out by comparing the 
right hand sides of the above ODEs. 
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